! include "blas.f90"

module PSTDDM
	implicit none
contains

	subroutine pstddm_gmres_conk_square(dpar,ipar,a_fun_full,c_fun_full, Fv, u, maxrelres, maxiter, relres, iter)
		! get the solution without zero boundary points of the PML equation: \divergenc(a_fun*\gradient(u)) + c_fun*u = Fv
		! by our PSTDDM as a preconditioner
		! dpar(1:2), the left buttom vertice of the whole domain
		! dpar(3), the width of PMLayer in x direction
		! dpar(4), the width of PMLayer in y direction
		! ipar(1), the number of intervals in dx_pml=dpar(3)
		! ipar(2), the number of intervals in dy_pml=dpar(4)
		! ipar(3), the multiple dx_pmls of the width of internal sub-domain in x direction
		! ipar(4), the multiple dy_pmls of the width of internal sub-domain in y direction
		! ipar(5), the number of layers in x direction
		! ipar(6), the number of layers in y direction
		! maxiter(1), the whole max iteration
		! maxiter(2), the max iteration in each restartion
		use sparseoperator
		implicit none
		interface
			subroutine a_fun_full(x,y,int_domain,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun_full
			subroutine c_fun_full(x,y,int_domain,c)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun_full
		end interface
		double precision, intent(in) :: dpar(4)
		integer, intent(in) :: ipar(6)
		complex(kind=8), dimension(:), intent(inout) :: Fv
		complex(kind=8), dimension(:), intent(inout) :: u
		double precision, intent(in) :: maxrelres
		integer, intent(in) :: maxiter(2)
		double precision, intent(inout) :: relres
		integer, intent(inout) :: iter

		double precision :: O(2), width, height, dx_pml, dy_pml, dx_seg, dy_seg
		integer :: m, n, m_pml, n_pml, m_pml_seg, n_pml_seg, Ny, Mx
		complex(kind=8), dimension(:), allocatable :: a
		integer, dimension(:), allocatable :: ia, ja
		double precision :: int_domain(4)
		integer :: dofs, nnz
		integer :: err


		O = dpar(1:2)
		dx_pml = dpar(3)
		dy_pml = dpar(4)
		m_pml = ipar(1)
		n_pml = ipar(2)
		m_pml_seg = ipar(3)
		n_pml_seg = ipar(4)
		Mx = ipar(5)
		Ny = ipar(6)

		width = 2.d0*dx_pml + Mx*m_pml_seg*dx_pml
		m = 2*m_pml + Mx*m_pml_seg*m_pml
		height = 2.d0*dy_pml + Ny*n_pml_seg*dy_pml
		n = 2*n_pml + Ny*n_pml_seg*n_pml

		call pstddm_getdn(m,n,dofs, nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(dofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		int_domain = (/ O(1)+dx_pml, O(1)+width-dx_pml, &
				&	O(2)+dy_pml, O(2)+height-dy_pml /)
		call pstddm_zstiffness(O,width,height,m,n, a_fun, c_fun, a, ia, ja)

		call sp_zgmres_precondition(maxiter, maxrelres, a,ia,ja,Fv, u, iter, relres, prefun)
		! 				call prefun(Fv,u)
		! 		call sp_zpardiso(13,0,a,ia,ja,Fv,u)


		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
	contains
		subroutine prefun(x,u)
			complex(kind=8), dimension(:), intent(in) :: x
			complex(kind=8), dimension(:), intent(inout) :: u
			call pstddm_conk_square(dpar,ipar,a_fun_full,c_fun_full,x,u)
		end subroutine prefun
		subroutine a_fun(x,y,a1,a2)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			call a_fun_full(x,y,int_domain,a1,a2)
		end subroutine a_fun
		subroutine c_fun(x,y,c)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: c
			call c_fun_full(x,y,int_domain,c)
		end subroutine c_fun
	end subroutine pstddm_gmres_conk_square

	subroutine pstddm_conk_square(dpar,ipar,a_fun_full,c_fun_full,Fv,u)
		! get the solution without zero boundary points of the PML equation: \divergenc(a_fun*\gradient(u)) + c_fun*u = Fv
		! by our PSTDDM as a direct solver
		! dpar(1:2), the left buttom vertice of the whole domain
		! dpar(3), the width of PMLayer in x direction
		! dpar(4), the width of PMLayer in y direction
		! ipar(1), the number of intervals in dx_pml=dpar(3)
		! ipar(2), the number of intervals in dy_pml=dpar(4)
		! ipar(3), the multiple dx_pmls of the width of internal sub-domain in x direction
		! ipar(4), the multiple dy_pmls of the width of internal sub-domain in y direction
		! ipar(5), the number of layers in x direction
		! ipar(6), the number of layers in y direction
		use sparseoperator
		implicit none
		interface
			subroutine a_fun_full(x,y,int_domain,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun_full
			subroutine c_fun_full(x,y,int_domain,c)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun_full
		end interface
		double precision, intent(in) :: dpar(4)
		integer, intent(in) :: ipar(5)
		complex(kind=8), dimension(:), intent(in) :: Fv
		complex(kind=8), dimension(:), intent(inout) :: u

		double precision :: O(2), width, height, dx_pml, dy_pml, dx_int,dx_seg, dy_seg, height_sub
		integer :: m_pml, n_pml, m_int, m_seg, n_seg, m, n, n_sub, m_pml_int, m_pml_seg, n_pml_seg, Mx, Ny
		integer :: err, ii, jj, ldofs, iny

		integer, dimension(:), allocatable :: ia, ja
		complex(kind=8), dimension(:), allocatable :: a, F_sub, u_sub

		double precision :: int_domain(4), hx, hy
		double precision, dimension(:), allocatable :: bt, dtmp

		integer :: mtype, msglvl, nnz

		double precision :: dparx(4)
		integer :: iparx(5)

		integer :: ios


		O = dpar(1:2)
		dx_pml = dpar(3)
		dy_pml = dpar(4)

		m_pml = ipar(1)
		n_pml = ipar(2)
		m_pml_seg = ipar(3)
		n_pml_seg = ipar(4)
		Mx = ipar(5)
		Ny = ipar(6)

		m_pml_int = Mx*m_pml_seg

		dx_int = m_pml_int * dx_pml
		width = 2.d0*dx_pml + dx_int
		dy_seg = n_pml_seg * dy_pml
		height = 2.d0*dy_pml + Ny*dy_seg
		height_sub = 2.d0*dy_pml + 2.d0*dy_seg

		m_int = m_pml_int * m_pml
		m = 2*m_pml + m_int
		n_seg = n_pml_seg * n_pml
		n_sub = 2*n_pml + 2*n_seg
		n = 2*n_pml + Ny*n_seg

		hx = width/m
		hy = height/n

		int_domain = (/ O(1)+dx_pml, O(1)+width-dx_pml, O(2)+dy_pml, O(2)+height_sub-dy_pml /)
		call pstddm_getdn(m,n_sub,ldofs,nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(ldofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		call pstddm_zstiffness(O,width,height_sub,m,n_sub, a_fun, c_fun, a, ia, ja)


		allocate(bt(ldofs), stat=err)
		if (err /= 0) print *, "bt: Allocation request denied"
		allocate(dtmp(ldofs), stat=err)
		if (err /= 0) print *, "dtmp: Allocation request denied"
		forall(ii=1:m-1,jj=1:n_sub-1)
			dtmp(ii+(jj-1)*(m-1)) = O(2) + jj*hy
		end forall
		call betai(dtmp, O(2)+dy_pml+dy_seg, O(2)+dy_pml+2.d0*dy_seg, bt, '+')


		allocate(F_sub(ldofs), stat=err)
		if (err /= 0) print *, "F_sub: Allocation request denied"
		allocate(u_sub(ldofs), stat=err)
		if (err /= 0) print *, "u_sub: Allocation request denied"
		F_sub = 0.d0
		u_sub = 0.d0
		u = 0.d0
		do iny = 1, Ny-1, 1
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				F_sub(ii+(jj-1)*(m-1)) = F_sub(ii+(jj-1)*(m-1)) &
					&	+ Fv( ii+(jj-1+(iny-1)*n_seg)*(m-1) )
			end forall
			if ( iny.eq.Ny-1 ) then
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1+n_seg)*(m-1)) = Fv( ii+(jj-1+iny*n_seg)*(m-1) )
				end forall
			end if
			dparx = (/ O(1), O(2)+(iny-1.d0)*dy_seg, dx_pml, dy_pml /)
			iparx = (/ m_pml, n_pml, m_pml_seg, 2*n_pml_seg, Mx /)
			call pstddm_conk_x(dparx,iparx,a_fun_full,c_fun_full,F_sub,u_sub)
			u_sub = -u_sub
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				u(ii+(jj-1+iny*n_seg)*(m-1)) = u_sub(ii+(jj-1+n_seg)*(m-1))
			end forall
			if ( iny.eq.Ny-1 ) then
				forall(ii=1:m-1, jj=n_pml+2*n_seg:n_sub-1 )
					u( ii + ( jj-1+ (Ny-2)*n_seg ) *(m-1) ) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
			if ( iny.lt.Ny-1 ) then
				F_sub = u_sub*bt
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1)*(m-1)) = u_sub(ii+(jj-1+n_seg)*(m-1))
				end forall
			end if
		end do

		call betai(dtmp, O(2)+dy_pml, O(2)+dy_pml+dy_seg, bt, '-')
		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		F_sub = 0.d0
		do iny = Ny, 2, -1
			forall(ii=1:m-1,jj=n_pml:n_pml+2*n_seg-1)
				F_sub(ii+(jj-1)*(m-1)) = F_sub(ii+(jj-1)*(m-1)) &
					& + Fv(ii+(jj-1+(iny-2)*n_seg)*(m-1))
			end forall
			dparx = (/ O(1), O(2)+(iny-2.d0)*dy_seg, dx_pml, dy_pml /)
			iparx = (/ m_pml, n_pml, m_pml_seg, 2*n_pml_seg, Mx /)
			call pstddm_conk_x(dparx,iparx,a_fun_full,c_fun_full,F_sub,u_sub)
			u_sub = - u_sub
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				u(ii+(jj-1+(iny-2)*n_seg)*(m-1)) = u(ii+(jj-1+(iny-2)*n_seg)*(m-1)) &
					&	+ u_sub(ii+(jj-1)*(m-1))
			end forall
			if ( iny.eq.2 ) then
				forall(ii=1:m-1,jj=1:n_pml+n_seg-1)
					u(ii+(jj-1)*(m-1)) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
			if ( iny.gt.2 ) then
				F_sub = bt*u_sub
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1+n_seg)*(m-1)) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
		end do

		u = -u;

		if (allocated(u_sub)) deallocate(u_sub, stat=err)
		if (err /= 0) print *, "u_sub: Deallocation request denied"

		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		if (allocated(bt)) deallocate(bt, stat=err)
		if (err /= 0) print *, "bt: Deallocation request denied"

		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
		if (allocated(F_sub)) deallocate(F_sub, stat=err)
		if (err /= 0) print *, "F_sub: Deallocation request denied"
	contains
		subroutine a_fun(x,y,a1,a2)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			call a_fun_full(x,y,int_domain,a1,a2)
		end subroutine a_fun

		subroutine c_fun(x,y,c)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: c
			call c_fun_full(x,y,int_domain,c)
		end subroutine c_fun

		subroutine betai(y,y1,y2,bt,drt)
			double precision, dimension(:), intent(in) :: y
			double precision, intent(in) :: y1, y2
			double precision, dimension(:), intent(inout) :: bt
			character(len=1), intent(in) :: drt

			double precision :: ly, dy

			integer :: ii, sy

			ly = 8.d0
			dy = y2-y1
			ly = dy / ly

			sy = size(y)
			if ( drt.eq.'+' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 1.d0
					elseif (y(ii).ge.y2-ly ) then
						bt(ii) = 0.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = ( bt(ii)-1.d0 )**2
					end if
				end do
			end if
			if ( drt.eq.'-' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 0.d0
					elseif ( y(ii).ge.y2-ly ) then
						bt(ii) = 1.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = bt(ii)*(2.d0-bt(ii))
					end if
				end do

			end if
		end subroutine betai
	end subroutine pstddm_conk_square

	subroutine pstddm_conk_x(dpar,ipar,a_fun_full,c_fun_full,Fv,u)
		! get the solution without zero boundary points of the PML equation: \divergenc(a_fun*\gradient(u)) + c_fun*u = Fv
		! by our PSTDDM as a direct solver
		! dpar(1:2), the left buttom vertice of the whole domain
		! dpar(3), the width of PMLayer in x direction
		! dpar(4), the width of PMLayer in y direction
		! ipar(1), the number of intervals in dx_pml=dpar(3)
		! ipar(2), the number of intervals in dy_pml=dpar(4)
		! ipar(3), the multiple dx_pmls of the width of internal sub-domain in x direction
		! ipar(4), the multiple dy_pmls of the width of internal domain in y direction
		! ipar(5), the number of layers in x direction
		use sparseoperator
		implicit none
		interface
			subroutine a_fun_full(x,y,int_domain,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun_full
			subroutine c_fun_full(x,y,int_domain,c)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun_full
		end interface
		double precision, intent(in) :: dpar(4)
		integer, intent(in) :: ipar(5)
		complex(kind=8), dimension(:), intent(in) :: Fv
		complex(kind=8), dimension(:), intent(inout) :: u

		double precision :: O(2), width, height, dx_pml, dy_pml, dx_seg, dy_int, width_sub
		integer :: m_pml, n_pml, m_seg, n_int, m, n, m_sub, m_pml_seg, n_pml_int, Mx
		integer :: err, ii, jj, ldofs, imx

		integer, dimension(:), allocatable :: ia, ja
		complex(kind=8), dimension(:), allocatable :: a, F_sub, u_sub

		double precision :: int_domain(4), hx, hy
		double precision, dimension(:), allocatable :: bt, dtmp

		integer :: mtype, msglvl, nnz


		O = dpar(1:2)
		dx_pml = dpar(3)
		dy_pml = dpar(4)

		m_pml = ipar(1)
		n_pml = ipar(2)
		m_pml_seg = ipar(3)
		n_pml_int = ipar(4)
		Mx = ipar(5)

		dx_seg = m_pml_seg * dx_pml
		width = 2.d0*dx_pml + Mx*dx_seg
		dy_int = n_pml_int * dy_pml
		height = 2.d0*dy_pml + dy_int
		width_sub = 2.d0*dx_pml + 2.d0*dx_seg

		m_seg = m_pml_seg * m_pml
		m = 2*m_pml + Mx*m_seg
		m_sub = 2*m_pml + 2*m_seg
		n_int = n_pml_int * n_pml

		n = 2*n_pml + n_int

		hx = width/m
		hy = height/n

		int_domain = (/ O(1)+dx_pml, O(1)+width_sub-dx_pml, O(2)+dy_pml, O(2)+height-dy_pml /)
		call pstddm_getdn(m_sub,n,ldofs,nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(ldofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		call pstddm_zstiffness(O,width_sub,height,m_sub,n, a_fun, c_fun, a, ia, ja)

		allocate(bt(ldofs), stat=err)
		if (err /= 0) print *, "bt: Allocation request denied"
		allocate(dtmp(ldofs), stat=err)
		if (err /= 0) print *, "dtmp: Allocation request denied"
		forall(ii=1:m_sub-1,jj=1:n-1)
			dtmp(ii+(jj-1)*(m_sub-1)) = O(1) + ii*hx
		end forall
		call betai(dtmp, O(1)+dx_pml+dx_seg, O(1)+dx_pml+2.d0*dx_seg, bt, '+')


		allocate(F_sub(ldofs), stat=err)
		if (err /= 0) print *, "F_sub: Allocation request denied"
		allocate(u_sub(ldofs), stat=err)
		if (err /= 0) print *, "u_sub: Allocation request denied"
		F_sub = 0.d0
		u_sub = 0.d0
		u = 0.d0
		do imx = 1, Mx-1, 1
			forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
				F_sub(ii+(jj-1)*(m_sub-1)) = F_sub(ii+(jj-1)*(m_sub-1)) &
					&	+ Fv( ii+(jj-1)*(m-1)+(imx-1)*m_seg )
			end forall
			if ( imx.eq.Mx-1 ) then
				forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
					F_sub(ii+(jj-1)*(m_sub-1)+m_seg) = Fv( ii+(jj-1)*(m-1)+imx*m_seg )
				end forall
			end if
			mtype = 13	! complex and nonsymmetric
			msglvl = 0	! do not print all the pardiso message
			call sp_zpardiso(mtype,msglvl,a,ia,ja,F_sub,u_sub)
			u_sub = -u_sub
			forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
				u(ii+(jj-1)*(m-1)+imx*m_seg) = u_sub(ii+(jj-1)*(m_sub-1)+m_seg)
			end forall
			if ( imx.eq.Mx-1 ) then
				forall(ii=m_pml+2*m_seg:m_sub-1, jj=1:n-1 )
					u( ii+(jj-1)*(m-1)+(Mx-2)*m_seg) = u_sub(ii+(jj-1)*(m_sub-1))
				end forall
			end if
			if ( imx.lt.Mx-1 ) then
				F_sub = u_sub*bt
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
					F_sub(ii+(jj-1)*(m_sub-1)) = u_sub(ii+(jj-1)*(m_sub-1)+m_seg)
				end forall
			end if
		end do

		call betai(dtmp, O(1)+dx_pml, O(1)+dx_pml+dx_seg, bt, '-')
		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		F_sub = 0.d0
		do imx = Mx, 2, -1
			forall(ii=m_pml:m_pml+2*m_seg-1,jj=1:n-1)
				F_sub(ii+(jj-1)*(m_sub-1)) = F_sub(ii+(jj-1)*(m_sub-1)) &
					& + Fv(ii+(jj-1)*(m-1)+(imx-2)*m_seg)
			end forall
			mtype = 13	! complex and nonsymmetric
			msglvl = 0	! do not print all the pardiso message
			call sp_zpardiso(mtype,msglvl,a,ia,ja,F_sub,u_sub)
			u_sub = - u_sub
			forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
				u(ii+(jj-1)*(m-1)+(imx-2)*m_seg) = u(ii+(jj-1)*(m-1)+(imx-2)*m_seg) &
					&	+ u_sub(ii+(jj-1)*(m_sub-1))
			end forall
			if ( imx.eq.2 ) then
				forall(ii=1:m_pml+m_seg-1,jj=1:n-1)
					u(ii+(jj-1)*(m-1)) = u_sub(ii+(jj-1)*(m_sub-1))
				end forall
			end if
			if ( imx.gt.2 ) then
				F_sub = bt*u_sub
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=m_pml:m_pml+m_seg-1, jj=1:n-1 )
					F_sub(ii+(jj-1)*(m_sub-1)+m_seg) = u_sub(ii+(jj-1)*(m_sub-1))
				end forall
			end if
		end do

		u = -u;

		if (allocated(u_sub)) deallocate(u_sub, stat=err)
		if (err /= 0) print *, "u_sub: Deallocation request denied"

		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		if (allocated(bt)) deallocate(bt, stat=err)
		if (err /= 0) print *, "bt: Deallocation request denied"

		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
		if (allocated(F_sub)) deallocate(F_sub, stat=err)
		if (err /= 0) print *, "F_sub: Deallocation request denied"
	contains
		subroutine a_fun(x,y,a1,a2)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			call a_fun_full(x,y,int_domain,a1,a2)
		end subroutine a_fun

		subroutine c_fun(x,y,c)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: c
			call c_fun_full(x,y,int_domain,c)
		end subroutine c_fun

		subroutine betai(y,y1,y2,bt,drt)
			double precision, dimension(:), intent(in) :: y
			double precision, intent(in) :: y1, y2
			double precision, dimension(:), intent(inout) :: bt
			character(len=1), intent(in) :: drt

			double precision :: ly, dy

			integer :: ii, sy

			ly = 8.d0
			dy = y2-y1
			ly = dy / ly

			sy = size(y)
			if ( drt.eq.'+' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 1.d0
					elseif (y(ii).ge.y2-ly ) then
						bt(ii) = 0.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = ( bt(ii)-1.d0 )**2
					end if
				end do
			end if
			if ( drt.eq.'-' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 0.d0
					elseif ( y(ii).ge.y2-ly ) then
						bt(ii) = 1.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = bt(ii)*(2.d0-bt(ii))
					end if
				end do

			end if
		end subroutine betai
	end subroutine pstddm_conk_x

	subroutine pstddm_gmres_conk_y(dpar,ipar,a_fun_full,c_fun_full, Fv, u, maxrelres, maxiter, relres, iter)
		! get the solution without zero boundary points of the PML equation: \divergenc(a_fun*\gradient(u)) + c_fun*u = Fv
		! by our PSTDDM as a preconditioner of GMRES
		! dpar(1:2), the left buttom vertice of the whole domain
		! dpar(3), the width of PMLayer in x direction
		! dpar(4), the width of PMLayer in y direction
		! ipar(1), the number of intervals in dx_pml=dpar(3)
		! ipar(2), the number of intervals in dy_pml=dpar(4)
		! ipar(3), the multiple dx_pmls of the width of internal domain in x direction
		! ipar(4), the multiple dy_pmls of the width of internal sub-domain in y direction
		! ipar(5), the number of layers in y direction
		! maxiter(1), the whole max iteration
		! maxiter(2), the max iteration in each restartion
		use sparseoperator
		implicit none
		interface
			subroutine a_fun_full(x,y,int_domain,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun_full
			subroutine c_fun_full(x,y,int_domain,c)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun_full
		end interface
		double precision, intent(in) :: dpar(4)
		integer, intent(in) :: ipar(5)
		complex(kind=8), dimension(:), intent(inout) :: Fv
		complex(kind=8), dimension(:), intent(inout) :: u
		double precision, intent(in) :: maxrelres
		integer, intent(in) :: maxiter(2)
		double precision, intent(inout) :: relres
		integer, intent(inout) :: iter

		double precision :: O(2), width, height, dx_pml, dy_pml, dx_int, dy_seg
		integer :: m, n, m_pml, n_pml, m_pml_int, n_pml_seg, Ny
		complex(kind=8), dimension(:), allocatable :: a
		integer, dimension(:), allocatable :: ia, ja
		double precision :: int_domain(4)
		integer :: dofs, nnz
		integer :: err


		O = dpar(1:2)
		dx_pml = dpar(3)
		dy_pml = dpar(4)
		m_pml = ipar(1)
		n_pml = ipar(2)
		m_pml_int = ipar(3)
		n_pml_seg = ipar(4)
		Ny = ipar(5)

		width = 2.d0*dx_pml + m_pml_int*dx_pml
		m = 2*m_pml + m_pml_int*m_pml
		height = 2.d0*dy_pml + Ny*n_pml_seg*dy_pml
		n = 2*n_pml + Ny*n_pml_seg*n_pml

		call pstddm_getdn(m,n,dofs, nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(dofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		int_domain = (/ O(1)+dx_pml, O(1)+width-dx_pml, &
				&	O(2)+dy_pml, O(2)+height-dy_pml /)
		call pstddm_zstiffness(O,width,height,m,n, a_fun, c_fun, a, ia, ja)

		call sp_zgmres_precondition(maxiter, maxrelres, a,ia,ja,Fv, u, iter, relres, prefun)

		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
	contains
		subroutine prefun(x,u)
			complex(kind=8), dimension(:), intent(in) :: x
			complex(kind=8), dimension(:), intent(inout) :: u
			call pstddm_conk_y(dpar,ipar,a_fun_full,c_fun_full,x,u)
		end subroutine prefun
		subroutine a_fun(x,y,a1,a2)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			call a_fun_full(x,y,int_domain,a1,a2)
		end subroutine a_fun
		subroutine c_fun(x,y,c)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: c
			call c_fun_full(x,y,int_domain,c)
		end subroutine c_fun
	end subroutine pstddm_gmres_conk_y

	subroutine pstddm_conk_y(dpar,ipar,a_fun_full,c_fun_full,Fv,u)
		! get the solution without zero boundary points of the PML equation: \divergenc(a_fun*\gradient(u)) + c_fun*u = Fv
		! by our PSTDDM as a direct solver
		! dpar(1:2), the left buttom vertice of the whole domain
		! dpar(3), the width of PMLayer in x direction
		! dpar(4), the width of PMLayer in y direction
		! ipar(1), the number of intervals in dx_pml=dpar(3)
		! ipar(2), the number of intervals in dy_pml=dpar(4)
		! ipar(3), the multiple dx_pmls of the width of internal domain in x direction
		! ipar(4), the multiple dy_pmls of the width of internal sub-domain in y direction
		! ipar(5), the number of layers in y direction
		use sparseoperator
		implicit none
		interface
			subroutine a_fun_full(x,y,int_domain,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun_full
			subroutine c_fun_full(x,y,int_domain,c)
				double precision, dimension(:,:), intent(in) :: x, y
				double precision, intent(in) :: int_domain(4)
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun_full
		end interface
		double precision, intent(in) :: dpar(4)
		integer, intent(in) :: ipar(5)
		complex(kind=8), dimension(:), intent(in) :: Fv
		complex(kind=8), dimension(:), intent(inout) :: u

		double precision :: O(2), width, height, dx_pml, dy_pml, dx_int, dy_seg, height_sub
		integer :: m_pml, n_pml, m_int, n_seg, m, n, n_sub, m_pml_int, n_pml_seg, Ny
		integer :: err, ii, jj, ldofs, iny

		integer, dimension(:), allocatable :: ia, ja
		complex(kind=8), dimension(:), allocatable :: a, F_sub, u_sub

		double precision :: int_domain(4), hx, hy
		double precision, dimension(:), allocatable :: bt, dtmp

		integer :: mtype, msglvl, nnz
		integer :: ios


		O = dpar(1:2)
		dx_pml = dpar(3)
		dy_pml = dpar(4)

		m_pml = ipar(1)
		n_pml = ipar(2)
		m_pml_int = ipar(3)
		n_pml_seg = ipar(4)
		Ny = ipar(5)

		dx_int = m_pml_int * dx_pml
		width = 2.d0*dx_pml + dx_int
		dy_seg = n_pml_seg * dy_pml
		height = 2.d0*dy_pml + Ny*dy_seg
		height_sub = 2.d0*dy_pml + 2.d0*dy_seg

		m_int = m_pml_int * m_pml
		m = 2*m_pml + m_int
		n_seg = n_pml_seg * n_pml
		n_sub = 2*n_pml + 2*n_seg
		n = 2*n_pml + Ny*n_seg

		hx = width/m
		hy = height/n

		int_domain = (/ O(1)+dx_pml, O(1)+width-dx_pml, O(2)+dy_pml, O(2)+height_sub-dy_pml /)
		call pstddm_getdn(m,n_sub,ldofs,nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(ldofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		call pstddm_zstiffness(O,width,height_sub,m,n_sub, a_fun, c_fun, a, ia, ja)

		allocate(bt(ldofs), stat=err)
		if (err /= 0) print *, "bt: Allocation request denied"
		allocate(dtmp(ldofs), stat=err)
		if (err /= 0) print *, "dtmp: Allocation request denied"
		forall(ii=1:m-1,jj=1:n_sub-1)
			dtmp(ii+(jj-1)*(m-1)) = O(2) + jj*hy
		end forall
		call betai(dtmp, O(2)+dy_pml+dy_seg, O(2)+dy_pml+2.d0*dy_seg, bt, '+')


		allocate(F_sub(ldofs), stat=err)
		if (err /= 0) print *, "F_sub: Allocation request denied"
		allocate(u_sub(ldofs), stat=err)
		if (err /= 0) print *, "u_sub: Allocation request denied"
		F_sub = 0.d0
		u_sub = 0.d0
		u = 0.d0

		ios = 0
		open(unit=1000, file='uf.txt', iostat=ios)
		if ( ios /= 0 ) stop "Error opening file name"

		do ii = 1, size(a), 1
			write(unit=1000, fmt='(2E20.12)') dreal(a(ii)), dimag(a(ii))
		end do

		do iny = 1, Ny-1, 1
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				F_sub(ii+(jj-1)*(m-1)) = F_sub(ii+(jj-1)*(m-1)) &
					&	+ Fv( ii+(jj-1+(iny-1)*n_seg)*(m-1) )
			end forall
			if ( iny.eq.Ny-1 ) then
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1+n_seg)*(m-1)) = Fv( ii+(jj-1+iny*n_seg)*(m-1) )
				end forall
			end if
			mtype = 13	! complex and nonsymmetric
			msglvl = 0	! do not print all the pardiso message
			call sp_zpardiso(mtype,msglvl,a,ia,ja,F_sub,u_sub)
			u_sub = -u_sub
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				u(ii+(jj-1+iny*n_seg)*(m-1)) = u_sub(ii+(jj-1+n_seg)*(m-1))
			end forall
			if ( iny.eq.Ny-1 ) then
				forall(ii=1:m-1, jj=n_pml+2*n_seg:n_sub-1 )
					u( ii + ( jj-1+ (Ny-2)*n_seg ) *(m-1) ) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
			if ( iny.lt.Ny-1 ) then
				F_sub = u_sub*bt
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1)*(m-1)) = u_sub(ii+(jj-1+n_seg)*(m-1))
				end forall
			end if
		end do

		call betai(dtmp, O(2)+dy_pml, O(2)+dy_pml+dy_seg, bt, '-')
		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		F_sub = 0.d0
		do iny = Ny, 2, -1
			forall(ii=1:m-1,jj=n_pml:n_pml+2*n_seg-1)
				F_sub(ii+(jj-1)*(m-1)) = F_sub(ii+(jj-1)*(m-1)) &
					& + Fv(ii+(jj-1+(iny-2)*n_seg)*(m-1))
			end forall
			mtype = 13	! complex and nonsymmetric
			msglvl = 0	! do not print all the pardiso message
			call sp_zpardiso(mtype,msglvl,a,ia,ja,F_sub,u_sub)
			u_sub = - u_sub
			forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
				u(ii+(jj-1+(iny-2)*n_seg)*(m-1)) = u(ii+(jj-1+(iny-2)*n_seg)*(m-1)) &
					&	+ u_sub(ii+(jj-1)*(m-1))
			end forall
			if ( iny.eq.2 ) then
				forall(ii=1:m-1,jj=1:n_pml+n_seg-1)
					u(ii+(jj-1)*(m-1)) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
			if ( iny.gt.2 ) then
				F_sub = bt*u_sub
				call mkl_zcsrgemv('n',ldofs,a,ia,ja,F_sub,u_sub)
				F_sub = 0.d0
				forall(ii=1:m-1, jj=n_pml:n_pml+n_seg-1 )
					F_sub(ii+(jj-1+n_seg)*(m-1)) = u_sub(ii+(jj-1)*(m-1))
				end forall
			end if
		end do

		u = -u;

		if (allocated(u_sub)) deallocate(u_sub, stat=err)
		if (err /= 0) print *, "u_sub: Deallocation request denied"

		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		if (allocated(bt)) deallocate(bt, stat=err)
		if (err /= 0) print *, "bt: Deallocation request denied"

		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
		if (allocated(F_sub)) deallocate(F_sub, stat=err)
		if (err /= 0) print *, "F_sub: Deallocation request denied"
	contains
		subroutine a_fun(x,y,a1,a2)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			call a_fun_full(x,y,int_domain,a1,a2)
		end subroutine a_fun

		subroutine c_fun(x,y,c)
			double precision, dimension(:,:), intent(in) :: x, y
			complex(kind=8), dimension(:,:), intent(inout) :: c
			call c_fun_full(x,y,int_domain,c)
		end subroutine c_fun

		subroutine betai(y,y1,y2,bt,drt)
			double precision, dimension(:), intent(in) :: y
			double precision, intent(in) :: y1, y2
			double precision, dimension(:), intent(inout) :: bt
			character(len=1), intent(in) :: drt

			double precision :: ly, dy

			integer :: ii, sy

			ly = 8.d0
			dy = y2-y1
			ly = dy / ly

			sy = size(y)
			if ( drt.eq.'+' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 1.d0
					elseif (y(ii).ge.y2-ly ) then
						bt(ii) = 0.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = ( bt(ii)-1.d0 )**2
					end if
				end do
			end if
			if ( drt.eq.'-' ) then
				do ii = 1, sy, 1
					if ( y(ii).le.y1+ly ) then
						bt(ii) = 0.d0
					elseif ( y(ii).ge.y2-ly ) then
						bt(ii) = 1.d0
					else
						bt(ii) = ((y(ii)-(y1+ly))/(dy-2.d0*ly))**2
						bt(ii) = bt(ii)*(2.d0-bt(ii))
					end if
				end do

			end if
		end subroutine betai
	end subroutine pstddm_conk_y

	subroutine pstddm_errors(O,width,height,m,n,u_fun,uxy_fun,u,m1,m2,n1,n2,errors)
		! get the errors in elements [m1:m2,n1:n2]
		! errors(1): the l2 norm of u_fun;
		! errors(2): the l2 error between u_fun and u
		! errors(3): the H1 norm of u_fun;
		! errors(4): the H1 error between u_fun and u
! 		use blas95
		implicit none
		interface
			subroutine u_fun(x,y,u)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: u
			end subroutine
			subroutine uxy_fun(x,y,ux,uy)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: ux, uy
			end subroutine
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n, m1, m2, n1, n2
		complex(kind=8), dimension(:), intent(in) :: u
		double precision, intent(inout) :: errors(4)


		double precision :: hx, hy, glpx(9), glpy(9), glw(9)
		double precision, dimension(:,:), allocatable :: x, y
		complex(kind=8), dimension(:,:), allocatable :: ue
		double precision, dimension(:), allocatable :: dtmp
		integer :: mr,nr, err, ii,jj
		double precision :: eps
		complex(kind=8) :: ddum(1,1)

		eps = 1.d-6

		hx = width/m
		hy = height/n

		glpx = (/	0.d0,	0.5d0,	1.d0,	&
			&	0.d0,	0.5d0,	1.d0,	&
			&	0.d0,	0.5d0,	1.d0 /)
		glpy = (/ 0.d0,	0.d0,	0.d0,	&
			&	0.5d0,	0.5d0,	0.5d0,	&
			&	1.d0,	1.d0,	1.d0 /)
		glw = (/	1.d0/36.d0,	1.d0/9.d0,	1.d0/36.d0,	&
			&	1.d0/9.d0,	4.d0/9.d0,	1.d0/9.d0,	&
			&	1.d0/36.d0,	1.d0/9.d0,	1.d0/36.d0 /) ! (/ 1.d0/6.d0, 2.d0/3.d0, 1.d0/6.d0 /)
		glw = glw * (hx*hy)
		mr = m2-m1+1
		nr = n2-n1+1

		if ( errors(1).gt.eps ) then
			allocate(x(mr*nr,9), stat=err)
			if (err /= 0) print *, "x: Allocation request denied"
			allocate(y(mr*nr,9), stat=err)
			if (err /= 0) print *, "y: Allocation request denied"

			forall(ii=1:mr,jj=1:nr)
				x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx*hx
				y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy*hy
			end forall

			allocate(ue(mr*nr,9), stat=err)
			if (err /= 0) print *, "ue: Allocation request denied"
			call u_fun(x,y,ue)


			allocate(dtmp(mr*nr), stat=err)
			if (err /= 0) print *, "dtmp: Allocation request denied"
			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw)
			end forall
			errors(1) = sqrt(sum(dtmp))
		end if

		if ( errors(2).gt.eps ) then
			if ( errors(1).lt.eps ) then
				allocate(x(mr*nr,9), stat=err)
				if (err /= 0) print *, "x: Allocation request denied"
				allocate(y(mr*nr,9), stat=err)
				if (err /= 0) print *, "y: Allocation request denied"

				forall(ii=1:mr,jj=1:nr)
					x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx*hx
					y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy*hy
				end forall

				allocate(ue(mr*nr,9), stat=err)
				if (err /= 0) print *, "ue: Allocation request denied"
				call u_fun(x,y,ue)
				allocate(dtmp(mr*nr), stat=err)
				if (err /= 0) print *, "dtmp: Allocation request denied"
			end if
			forall(ii=1:mr,jj=1:nr)
				ue(ii+(jj-1)*mr,:) = ue(ii+(jj-1)*mr,:)	&
					&	- u(m1-1+ii-1+(n1-1+jj-2)*(m-1))* (1.d0-glpx)*(1.d0-glpy)	&
					&	- u(m1-1+ii+(n1-1+jj-2)*(m-1))* glpx*(1.d0-glpy)	&
					&	- u(m1-1+ii-1+(n1-1+jj-1)*(m-1))* (1.d0-glpx)*glpy 	&
					&	- u(m1-1+ii+(n1-1+jj-1)*(m-1))* glpx*glpy
			end forall
			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw)
			end forall
			errors(2) = sqrt(sum(dtmp))
		end if

		if ( errors(3).gt.eps ) then
			if ( errors(1).lt.eps .and. errors(2).lt.eps ) then
				allocate(x(mr*nr,9), stat=err)
				if (err /= 0) print *, "x: Allocation request denied"
				allocate(y(mr*nr,9), stat=err)
				if (err /= 0) print *, "y: Allocation request denied"

				forall(ii=1:mr,jj=1:nr)
					x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx*hx
					y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy*hy
				end forall

				allocate(ue(mr*nr,9), stat=err)
				if (err /= 0) print *, "ue: Allocation request denied"
				allocate(dtmp(mr*nr), stat=err)
				if (err /= 0) print *, "dtmp: Allocation request denied"
			end if

			call uxy_fun(x,y,ue,ddum)

			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw)
			end forall
			errors(3) = sum(dtmp)
		end if

		if ( errors(4).gt.eps ) then
			if ( errors(1).lt.eps .and. errors(2).lt.eps .and. errors(3).lt.eps ) then
				allocate(x(mr*nr,9), stat=err)
				if (err /= 0) print *, "x: Allocation request denied"
				allocate(y(mr*nr,9), stat=err)
				if (err /= 0) print *, "y: Allocation request denied"

				forall(ii=1:mr,jj=1:nr)
					x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx*hx
					y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy*hy
				end forall
				allocate(ue(mr*nr,9), stat=err)
				if (err /= 0) print *, "ue: Allocation request denied"
				allocate(dtmp(mr*nr), stat=err)
				if (err /= 0) print *, "dtmp: Allocation request denied"
			end if
			if ( errors(3).lt.eps )	call uxy_fun(x,y,ue,ddum)

			forall(ii=1:mr,jj=1:nr)
				ue(ii+(jj-1)*mr,:) = ue(ii+(jj-1)*mr,:)	&
					&	- u(m1-1+ii-1+(n1-1+jj-2)*(m-1))* (-1.d0/hx)*(1.d0-glpy)	&
					&	- u(m1-1+ii+(n1-1+jj-2)*(m-1))* (1.d0/hx)*(1.d0-glpy)	&
					&	- u(m1-1+ii-1+(n1-1+jj-1)*(m-1))* (-1.d0/hx)*glpy 	&
					&	- u(m1-1+ii+(n1-1+jj-1)*(m-1))* (1.d0/hx)*glpy
			end forall

			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw )
			end forall
			errors(4) = sum(dtmp)
		end if

		if ( errors(3).gt.eps ) then
			if ( all(errors.lt.eps) ) then
				allocate(x(mr*nr,9), stat=err)
				if (err /= 0) print *, "x: Allocation request denied"
				allocate(y(mr*nr,9), stat=err)
				if (err /= 0) print *, "y: Allocation request denied"
				glpx = glpx * hx
				glpy = glpy * hy
				forall(ii=1:mr,jj=1:nr)
					x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx
				y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy
				end forall
				glpx = glpx / hx
				glpy = glpy / hy

				allocate(ue(mr*nr,9), stat=err)
				if (err /= 0) print *, "ue: Allocation request denied"
				allocate(dtmp(mr*nr), stat=err)
				if (err /= 0) print *, "dtmp: Allocation request denied"
			end if
			call uxy_fun(x,y,ddum,ue)
			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw)
			end forall
			errors(3) = sqrt( errors(3) + sum(dtmp) )
		end if

		if ( errors(4).gt.0.d0 ) then
			if ( all(errors.lt.eps) ) then
				allocate(x(mr*nr,9), stat=err)
				if (err /= 0) print *, "x: Allocation request denied"
				allocate(y(mr*nr,9), stat=err)
				if (err /= 0) print *, "y: Allocation request denied"
				glpx = glpx * hx
				glpy = glpy * hy
				forall(ii=1:mr,jj=1:nr)
					x(ii+(jj-1)*mr,:) = (O(1)+(m1-1+ii-1)*hx) + glpx
				y(ii+(jj-1)*mr,:) = (O(2)+(n1-1+jj-1)*hy) + glpy
				end forall
				glpx = glpx / hx
				glpy = glpy / hy
				allocate(ue(mr*nr,9), stat=err)
				if (err /= 0) print *, "ue: Allocation request denied"
				allocate(dtmp(mr*nr), stat=err)
				if (err /= 0) print *, "dtmp: Allocation request denied"
			end if
			if ( errors(3).lt.eps )	call uxy_fun(x,y,ddum,ue)
			forall(ii=1:mr,jj=1:nr)
				ue(ii+(jj-1)*mr,:) = ue(ii+(jj-1)*mr,:)	&
					&	- u(m1-1+ii-1+(n1-1+jj-2)*(m-1))* (1.d0-glpx)*(-1.d0/hy)	&
					&	- u(m1-1+ii+(n1-1+jj-2)*(m-1))* glpx*(-1.d0/hy)	&
					&	- u(m1-1+ii-1+(n1-1+jj-1)*(m-1))* (1.d0-glpx)/hy 	&
					&	- u(m1-1+ii+(n1-1+jj-1)*(m-1))* glpx/hy
			end forall
			forall(ii=1:mr,jj=1:nr)
				dtmp(ii+(jj-1)*mr) = dot_product( abs(ue(ii+(jj-1)*mr,:))**2, glw)
			end forall
			errors(4) = sqrt( errors(4) + sum(dtmp) )
		end if


		if (allocated(dtmp)) deallocate(dtmp, stat=err)
		if (err /= 0) print *, "dtmp: Deallocation request denied"

		if (allocated(ue)) deallocate(ue, stat=err)
		if (err /= 0) print *, "ue: Deallocation request denied"

		if (allocated(y)) deallocate(y, stat=err)
		if (err /= 0) print *, "y: Deallocation request denied"

		if (allocated(x)) deallocate(x, stat=err)
		if (err /= 0) print *, "x: Deallocation request denied"
	end subroutine pstddm_errors


	subroutine pstddm_zfdm(O,width,height,m,n, a_fun, c_fun, f_fun, u)
		! the numerical solution without zeros boundary of the finite difference method
		! for the problem:
		!			-\divergence(a_fun*\grandient(u)) + c_fun*u = f_fun
		! O, vector with size 2, the left bottum vertice
		! m, the intervals between (/ O(1), O(1)+width /)
		! n, the intervals between (/ O(2), O(2)+height /)
		use sparseoperator
		implicit none
		interface
			subroutine a_fun(x,y,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun
			subroutine c_fun(x,y,c)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun
			subroutine f_fun(x,y,f)
				double precision, dimension(:), intent(in) :: x, y
				complex(kind=8), dimension(:), intent(inout) :: f
			end subroutine
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		complex(kind=8), dimension(:), intent(inout) :: u

		complex(kind=8), dimension(:), allocatable :: a, Fv
		integer, dimension(:), allocatable :: ia, ja
		integer :: dofs, mtype, msglvl, err, nnz

		call pstddm_getdn(m,n,dofs, nnz)
		allocate(a(nnz), stat=err)
		if (err /= 0) print *, "a: Allocation request denied"
		allocate(ja(nnz), stat=err)
		if (err /= 0) print *, "ja: Allocation request denied"
		allocate(ia(dofs+1), stat=err)
		if (err /= 0) print *, "ia: Allocation request denied"
		call pstddm_zstiffness(O,width,height,m,n, a_fun, c_fun, a, ia, ja)

		allocate(Fv(dofs), stat=err)
		if (err /= 0) print *, "Fv: Allocation request denied"
		call pstddm_zrightvector(O,width,height,m,n, f_fun, Fv)

		mtype = 13	! complex and nonsymmetric
		msglvl = 0	! do not print all the pardiso message

		call sp_zpardiso(mtype,msglvl,a,ia,ja,Fv,u)

		if (allocated(a)) deallocate(a, stat=err)
		if (err /= 0) print *, "a: Deallocation request denied"
		if (allocated(ia)) deallocate(ia, stat=err)
		if (err /= 0) print *, "ia: Deallocation request denied"
		if (allocated(ja)) deallocate(ja, stat=err)
		if (err /= 0) print *, "ja: Deallocation request denied"
		if (allocated(Fv)) deallocate(Fv, stat=err)
		if (err /= 0) print *, "Fv: Deallocation request denied"
	end subroutine pstddm_zfdm

	subroutine pstddm_zrightvector(O,width,height,m,n, f_fun, Fv)
		! the right vector without zeros boundary of the finite difference method
		! O, vector with size 2, the left bottum vertice
		! m, the intervals between (/ O(1), O(1)+width /)
		! n, the intervals between (/ O(2), O(2)+height /)
		interface
			subroutine f_fun(x,y,f)
				double precision, dimension(:), intent(in) :: x, y
				complex(kind=8), dimension(:), intent(inout) :: f
			end subroutine
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		complex(kind=8), dimension(:), intent(inout) :: Fv

		double precision, dimension(:), allocatable :: x, y
		integer :: ii, jj, dofs, err
		double precision :: hx, hy

		hx = width/m
		hy = height/n
		dofs = (m-1)*(n-1)

		allocate(x(dofs), stat=err)
		if (err /= 0) print *, "x: Allocation request denied"
		allocate(y(dofs), stat=err)
		if (err /= 0) print *, "y: Allocation request denied"
		forall(ii=1:m-1,jj=1:n-1)
			x(ii+(jj-1)*(m-1)) = O(1)+ii*hx
			y(ii+(jj-1)*(m-1)) = O(2)+jj*hy
		end forall

		call f_fun(x,y,Fv)

		if (allocated(y)) deallocate(y, stat=err)
		if (err /= 0) print *, "y: Deallocation request denied"
		if (allocated(x)) deallocate(x, stat=err)
		if (err /= 0) print *, "x: Deallocation request denied"
	end subroutine pstddm_zrightvector


	subroutine pstddm_zstiffness(O,width,height,m,n, a_fun, c_fun, a, ia, ja)
		! the stiffness matrix without zeros boundary of the finite difference method
		! for the problem:
		!			-\nabla\cdot(a_fun*\nabla u) + c_fun*u = f_fun
		! O, vector with size 2, the left bottum vertice
		! m, the intervals between (/ O(1), O(1)+width /)
		! n, the intervals between (/ O(2), O(2)+height /)
		! a,ia,ja, the stiffness matrix in CSR format
		interface
			subroutine a_fun(x,y,a1,a2)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: a1, a2
			end subroutine a_fun
			subroutine c_fun(x,y,c)
				double precision, dimension(:,:), intent(in) :: x, y
				complex(kind=8), dimension(:,:), intent(inout) :: c
			end subroutine c_fun
		end interface
		double precision, intent(in) :: O(2), width, height
		integer, intent(in) :: m, n
		complex(kind=8), dimension(:), intent(inout) :: a
		integer, dimension(:), intent(inout) :: ia, ja

		double precision, dimension(:,:), allocatable :: x, y
		complex(kind=8), dimension(:,:), allocatable :: a1, a2, c
		integer :: ii, jj, dofs, nnz, err
		double precision :: hx, hy

		hx = width/m
		hy = height/n


		allocate(x(m+1,n+1), stat=err)
		if (err /= 0) print *, "x: Allocation request denied"
		allocate(y(m+1,n+1), stat=err)
		if (err /= 0) print *, "y: Allocation request denied"
		forall(ii=1:m+1,jj=1:n+1)
			x(ii,jj) = O(1)+(ii-1)*hx
			y(ii,jj) = O(2)+(jj-1)*hy
		end forall
		allocate(a1(m+1,n+1), stat=err)
		if (err /= 0) print *, "a1: Allocation request denied"
		allocate(a2(m+1,n+1), stat=err)
		if (err /= 0) print *, "a2: Allocation request denied"
		allocate(c(m+1,n+1), stat=err)
		if (err /= 0) print *, "c: Allocation request denied"


		call a_fun(x,y,a1,a2)

		call c_fun(x,y,c)

		if (allocated(x)) deallocate(x, stat=err)
		if (err /= 0) print *, "x: Deallocation request denied"
		if (allocated(y)) deallocate(y, stat=err)
		if (err /= 0) print *, "y: Deallocation request denied"

		call pstddm_getdn(m,n,dofs, nnz)

		ia(1) = 1
		ia(2) = 3
		ia(m) = 3
		ia((n-2)*(m-1)+1+1) = 3
		ia(dofs+1) = 3
		forall(ii=2:m-2)
			ia(ii+1) = 4
			ia(ii+1+(n-2)*(m-1)) = 4
		end forall
		forall(jj=2:n-2)
			ia(2+(jj-1)*(m-1)) = 4
			ia(jj*(m-1)+1) = 4
		end forall
		forall(ii=2:m-2,jj=2:n-2)
			ia(ii+1+(jj-1)*(m-1)) = 5
		end forall
		do ii = 2, dofs+1, 1
			ia(ii) = ia(ii) + ia(ii-1)
		end do

		ja(ia(1):ia(2)-1) = (/ 1, 2, m /)
		ja(ia(m-1):ia(m)-1) = (/ m-2, m-1, 2*(m-1) /)
		ja(ia((n-2)*(m-1)+1):ia((n-2)*(m-1)+2)-1) = (/ (n-3)*(m-1)+1, (n-2)*(m-1)+1, (n-2)*(m-1)+2 /)
		ja(ia(dofs):ia(dofs+1)-1) = (/ dofs-(m-1), dofs-1, dofs /)
		forall(ii=2:m-2)
			ja(ia(ii):ia(ii+1)-1) = (/ ii-1, ii, ii+1, ii+m-1 /)
			ja(ia(ii+(n-2)*(m-1)):ia(ii+1+(n-2)*(m-1))-1) = (/ ii+(n-3)*(m-1), ii-1+(n-2)*(m-1), ii+(n-2)*(m-1), ii+1+(n-2)*(m-1) /)
		end forall
		forall(jj=2:n-2)
			ja(ia(1+(jj-1)*(m-1)):ia(2+(jj-1)*(m-1))-1) = (/ 1+(jj-2)*(m-1), 1+(jj-1)*(m-1), 2+(jj-1)*(m-1), 1+jj*(m-1) /)
			ja(ia(jj*(m-1)):ia(jj*(m-1)+1)-1) = (/ (jj-1)*(m-1), jj*(m-1)-1, jj*(m-1), (jj+1)*(m-1) /)
		end forall
		forall(ii=2:m-2,jj=2:n-2)
			ja(ia(ii+(jj-1)*(m-1)):ia(ii+1+(jj-1)*(m-1))-1) = (/ ii+(jj-2)*(m-1), ii-1+(jj-1)*(m-1), ii+(jj-1)*(m-1), ii+1+(jj-1)*(m-1), ii+jj*(m-1) /)
		end forall

		a1 = a1 / hx**2
		a2 = a2 / hy**2
		a(ia(1):ia(2)-1) = (/ a1(2,2)*2.d0+a2(2,2)*2.d0+c(2,2),	&
			&	-a1(2,2)-(a1(3,2)-a1(1,2))/4.d0,	&
			&	-a2(2,2)-(a2(2,3)-a2(2,1))/4.d0 /)

		a(ia(m-1):ia(m)-1) = (/ -a1(m,2)+(a1(m+1,2)-a1(m-1,2))/4.d0,	&
			&	a1(m,2)*2.d0+a2(m,2)*2.d0+c(m,2),	&
			&	-a2(m,2)-(a2(m,3)-a2(m,1))/4.d0 /)
		a(ia((n-2)*(m-1)+1):ia((n-2)*(m-1)+2)-1) = (/ -a2(2,n)+(a2(2,n+1)-a2(2,n-1))/4.d0,	&
			&	a1(2,n)*2.d0 + a2(2,n)*2.d0 + c(2,n), &
			&	-a1(2,n) - (a1(3,n)-a1(1,n))/4.d0 /)
		a(ia(dofs):ia(dofs+1)-1) = (/ -a2(m,n)+(a2(m,n+1)-a2(m,n-1))/4.d0,	&
			&	-a1(m,n)+(a1(m+1,n)-a1(m-1,n))/4.d0,	&
			&	a1(m,n)*2.d0+a2(m,n)*2.d0+c(m,n) /)
		forall(ii=2:m-2)
			a(ia(ii):ia(ii+1)-1) = (/ -a1(ii+1,2)+(a1(ii+2,2)-a1(ii,2))/4.d0,	&
				&	a1(ii+1,2)*2.d0+a2(ii+1,2)*2.d0+c(ii+1,2),	&
				&	-a1(ii+1,2)-(a1(ii+2,2)-a1(ii,2))/4.d0,	&
				&	-a2(ii+1,2)-(a2(ii+1,3)-a2(ii+1,1))/4.d0 /)
			a(ia(ii+(n-2)*(m-1)):ia(ii+1+(n-2)*(m-1))-1) = (/ -a2(ii+1,n)+(a2(ii+1,n+1)-a2(ii+1,n-1))/4.d0,	&
				&	-a1(ii+1,n)+(a1(ii+2,n)-a1(ii,n))/4.d0,	&
				&	a1(ii+1,n)*2.d0+a2(ii+1,n)*2.d0+c(ii+1,n),	&
				&	-a1(ii+1,n)-(a1(ii+2,n)-a1(ii,n))/4.d0 /)
		end forall
		forall(jj=2:n-2)
			a(ia(1+(jj-1)*(m-1)):ia(2+(jj-1)*(m-1))-1) = (/ -a2(2,jj+1)+(a2(2,jj+2)-a2(2,jj))/4.d0,	&
				&	a1(2,jj+1)*2.d0+a2(2,jj+1)*2.d0+c(2,jj+1),	&
				&	-a1(2,jj+1)-(a1(3,jj+1)-a1(1,jj+1))/4.d0,	&
				&	-a2(2,jj+1)-(a2(2,jj+2)-a2(2,jj))/4.d0 /)
			a(ia(jj*(m-1)):ia(jj*(m-1)+1)-1) = (/ -a2(m,jj+1)+(a2(m,jj+2)-a2(m,jj))/4.d0,	&
				&	-a1(m,jj+1)+(a1(m+1,jj+1)-a1(m-1,jj+1))/4.d0,	&
				&	a1(m,jj+1)*2.d0+a2(m,jj+1)*2.d0+c(m,jj+1),	&
				&	-a2(m,jj+1)-(a2(m,jj+2)-a2(m,jj))/4.d0 /)
		end forall
		forall(ii=2:m-2,jj=2:n-2)
			a(ia(ii+(jj-1)*(m-1)):ia(ii+1+(jj-1)*(m-1))-1) = (/ -a2(ii+1,jj+1)+(a2(ii+1,jj+2)-a2(ii+1,jj))/4.d0,	&
				&	-a1(ii+1,jj+1)+(a1(ii+2,jj+1)-a1(ii,jj+1))/4.d0,	&
				&	a1(ii+1,jj+1)*2.d0+a2(ii+1,jj+1)*2.d0+c(ii+1,jj+1),	&
				&	-a1(ii+1,jj+1)-(a1(ii+2,jj+1)-a1(ii,jj+1))/4.d0,	&
				&	-a2(ii+1,jj+1)-(a2(ii+1,jj+2)-a2(ii+1,jj))/4.d0 /)
		end forall

		if ( allocated(x) ) deallocate(x)
		if ( allocated(y) ) deallocate(y)
		if ( allocated(a1) ) deallocate(a1)
		if ( allocated(a2) ) deallocate(a2)
		if ( allocated(c) ) deallocate(c)
	end subroutine pstddm_zstiffness

	subroutine pstddm_getdn(m,n,dofs, nnz)
		! get the numbers of the degree of freedoms and the non-zeros elements without zero boundary condition
		integer, intent(in) :: m, n
		integer, intent(inout) :: dofs
		integer,intent(inout), optional :: nnz
		if ( m.lt.3 .or. n.lt.3 ) then
			print *, "pstddm_getdn: m, n should be greater than or equal to 3"
		end if
		dofs = (m-1)*(n-1)
		if ( present(nnz) ) nnz = 3*4 + 4*2*(m-3+n-3) + 5*(m-3)*(n-3)
	end subroutine pstddm_getdn

end module PSTDDM